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Tight-binding (TB) molecular dynamics (MD) has emerged as a powerful method for investigating 
the atomic-scale structure of materials — in particular the interplay between structural and elec- 
tronic properties — bridging the gap between empirical methods which, while fast and efficient, lack 
transferability, and ab initio approaches which, because of excessive computational workload, suffer 
from limitations in size and run times. In this short review article, we examine several recent ap- 
plications of TBMD in the area of defects in covalently-bonded semiconductors and the amorphous 
phases of these materials. 
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I. INTRODUCTION 

As one will be able to judge from this special issue 
of Computational Materials Science, tight-binding (TB) 
molecular dynamics (MD) has evolved into a power- 
ful method for understanding material properties at the 
atomic level, offering a good compromise between empir- 
ical [D and first-principles approaches for describing 
the interactions between atoms. Indeed, empirical meth- 
ods lack transferability — model potentials are usually 
fitted to specific material properties in specific configura- 
tions, and often fail to properly describe situations and 
properties other than those on which they were fitted. 
In contrast, first-principles methods are transferable, but 
their computational workload is so great that only small 
systems (less than a hundred particles or so) can be dealt 
with on relatively short timescales (at most a hundred ps 
or so) on the fastest computers. 

With TBMD, and in particular the novel O(N) meth- 
ods (for a review, see, for instance, Ordejon's paper in 
this journal and Ref. it is possible to simulate sys- 
tems containing several hundred particles for a good frac- 
tion of a ns. This allows a number of interesting prob- 
lems to be addressed, as we illustrate below. In fact, the 
accuracy and power of the method can be enhanced con- 
siderably by combining it with other approaches, either 
empirical or first principles; we also give examples of this 
below. 

The purpose of the present review is to illustrate the 
scope of application of the TBMD method by means of 
selected examples. TB is a field that has been active for 
some time in the world of electronic structure calcula- 
tions Q , but only in recent years has it been coupled to 
MD, making it possible to study the interplay between 
structure and physical properties of materials. Thus it 



is possible, with TBMD, to investigate dynamical prop- 
erties per se, including relaxation, as well to optimise 
structural models, after which the electronic and other 
properties can be determined. 

We focus here on covalently-bonded semiconductors, 
mostly Si and the III-V's; carbon is the object of an- 
other article in this Journal. This review is not meant to 
be exhaustive — but rather illustrative — and thus nec- 
essarily incomplete; we therefore apologize to all whose 
work is not covered or mentioned. Two classes of prob- 
lems are examined: first, defects — be they localised 
or extended — and, second, amorphous covalent semi- 
conductors. TBMD has allowed significant progress to 
be made in the study of defects: because they are of 
quantum-mechanical nature, TB potentials are more re- 
liable and more accurate than empirical ones; at the same 
time, because they are semi-empirical, they allow larger 
systems to be studied on longer timescales than fully ab 
initio approaches. The same same applies to the study 
of amorphous materials, where the principal difficulty is 
in the proper description of the wide spectrum of highly- 
strained environments that are found in these materials. 
Of particular interest is the relation between structural 
and electronic properties which, evidently, cannot be de- 
rived from empirical models. Before discussing these ap- 
plications, we provide, for completeness, a short overview 
of the methodology; more details on TBMD can be found 
in, for instance, Ref. ||, as well as other articles in this 
issue; an excellent discussion of MD can be found in Ref. 
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II. TBMD 

TB is a standard method H for computing the elec- 
tronic properties of materials in terms of a set of pa- 
rameters describing the overlap between atomic orbitals 
on nearest-neighbour and, sometimes, second-nearest- 
neighbour sites. In order to carry out MD or static relax- 
ation calculations (i.e., to compute the forces), however, 
it is necessary to add to the total energy a repulsive term 
which includes a position-dependent electronic bonding 
energy as well as ionic contributions. The total energy 
for TBMD simulations can thus be written as 

Etot = 2_j + ^ C^rcp(|Rj — R-j I ) , 

i ij 

where the first term on the right is the "band-structure" 
energy, i.e., the quantum- mechanical bonding energy re- 
sulting from the overlap of atomic orbitals, and the sec- 
ond term is a two-body classical potential which ac- 
counts for all the other contributions to the total energy. 
The electronic eigenvalues are obtained from a simplified 
local-atomic-orbital representation of the Hamiltonian. 
The electronic wave-functions are generally expanded in 
terms of a reduced number of localised, orthonormal ba- 
sis functions 

A* 

where the coefficients are obtained by solving 

if 

In general, the basis set is restricted to the valence 
electron states. In the case of silicon, for example, one 
typically uses the single 3s and the three 3p orbitals — a 
much smaller basis set than would be needed in a plane- 
wave description of the electron wavefunctions. 

Different potentials differ in the way that the Hamilto- 
nian matrix elements H^ v are approximated, but also in 
the functional form of the two-body potential. The ma- 
trix elements — which depend a priori on distance and 
bond angles — are determined either from ab initio cal- 
culations or extracted from experiment. The most com- 
mon approximation consists in parametrising the over- 
lap integral in terms of a set of constants which decay as 
1 /r 2 , with a cutoff distance between the first and second- 
neighbour shells; a short cutoff distance speeds up the 
calculations but can be a source of problems in disor- 
dered materials where near-neighbour shells are not per- 
fectly separated. It is also possible to compute the matrix 
elements and the repulsive potential in a more accurate 
way using either the local-density approximation |J or 
the local-orbital representation of density-functional the- 
ory first introduced by Harris 0. These schemes, often 
dubbed "ab initio TB", constitute a trade-off between 
the accrued computational effort associated with them 



and the more accurate description of strained environ- 
ments which is of particular importance for disordered 
materials and defects. 



III. DEFECTS 

Defects play a major role in determining the physical 
properties of semiconductors H|. Even when present at 
low density, they affect deeply the electronic structure of 
the materials. This is true of point defects, but also of 
extended defects, especially in the case of quantum struc- 
tures. In spite of numerous experimental or theoretical 
studies, a complete picture of the structure of even the 
simplest defects (vacancies, interstitials, and small com- 
plexes of them, such as divacancies), in the most-studied 
semiconductor material — silicon — has not yet emerged. 
Defects, further, are not static objects, and can undergo 
diffusion at sufficiently high temperature. Again, little is 
known of such processes; the diffusion coefficient of H in 
Si, for instance, is still not understood in detail. 

TBMD calculations of defects have contributed to our 
understanding of their static and dynamic properties, but 
they have also been extremely useful in validating the 
models. Indeed, because they break the local symmetry, 
which often results in subtle relaxation and electronic ef- 
fects, defects are difficult to treat using empirical models 
and therefore serve as an excellent test of the ability of 
TB models to deal with low-symmetry situations. The 
details of the atomic structure of defects, however, are 
generally not known from experiment, and the tests must 
be against ab initio calculations, which themselves carry 
significant uncertainties because of limitations in size and 
computational load. 

Here we examine recents results on intrinsic (or native) 
defects. Extrinsic defects, i.e., impurities, are a difficult 
problem because of the additional complexity involved in 
constructing multi-species interactions. We nevertheless 
consider some exceptions, notably H diffusion in Si and 
B relaxation in Si. 



A. Intrinsic point defects in silicon 

1. Basics 

One of the first applications of TBMD was the study, 
by Wang et al. || , of native point defects in crystalline sil- 
icon using the TB model of Goodwin, Skinner and Petti- 
for (GSP) |l(| . Wang et al. have examined the formation 
energies of neutral monovacancies and self-interstitials as 
a function of the size of the simulation cell — up to 512 
atoms. The results are listed in Table |: One sees that 
size somewhat influences the formation energy, especially 
for the monovacancy, indicating that the distortion pat- 
tern of the defects cannot be accommodated fully by a 
small cell. The relaxed values of the formation energies 
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seem to fall close to results obtained in the local-density 
approximation (LDA) fl4|-|l6||, probably within the un- 
certainties inherent to both approaches. 

For the single vacancy, Wang et al. find a tetragonal 
Jahn- Teller distortion on top of the radial displacement 
of nearest-neighbour atoms, as was also observed by Song 
et al. jL7|. This agrees with the electron paramagnetic 
resonance measurements of Watkins 18 1, as well as the 
LDA calculations of Baraff et al. |L9| and, more recently, 
Seong and Lewis The LDA calculations, however, 

predict slightly larger total displacements than TB - 
0.4 versus about 0.3 A. The Jahn- Teller distortion is of 
quantum-mechanical origin and therefore not available 
from empirical models. 

Song et al. j!?) have used the GSP model to inves- 
tigate, in addition, the simple and split divacancies, as 
well as the Frenkel pair, all in their neutral charge states; 
the results as shown in Table |[ The formation energies 
are large and the relaxation energies — the difference 
in energy between unrclaxed and relaxed configurations, 
given in parenthesis in Table || — clearly non-negligible, 
typically representing a good fraction (30% or so) of the 
formation energy. Of course, this is accompanied by a sig- 
nificant change in volume during relaxation, and atomic 
displacements that can be as large as 1.25 A (in the case 
of the split divacancy |l7f| ). 

Within the GSP-TB model it is energetically 
favourable for two vacancies to "coalesce" , saving about 
1.68 eV in the process. Indeed, the formation energy of 
two isolated vacancies is 7.36 eV, dropping to 6.54 eV 
for the split divacancy, and to 5.68 for the simple diva- 
cancy. Thus, divacancies are expected to readily form 
and be relatively stable even at high temperatures. Like- 
wise, the formation of a Frenkel pair by a vacancy and 
an interstitial can reduce their total energy by as much 
as 1.43 eV - from 7.98 to 5.55 eV. 

The activation energy for diffusion is the sum of forma- 
tion and migration energies. The latter is the energy at 
the transition state between two equilibrium sites. Song 
et al. fT^I estimate the migration energy for the vacancy 
to be less than 1.0 eV so the activation energy, using the 
formation energy values discussed previously, must be 
less than 4.7 eV. Likewise, tetrahedral interstitials mi- 
grate via hexagonal sites with an energy of about 0.63 
eV, and thus the activation energy in this case would be 
of the order of 5.0 eV. 

As mentioned above, the atomic displacements for the 
monovacancy are predicted by the LDA to be slightly 
larger than the TB values. The opposite is true in the 
case of divacancies, where the LDA predicts displace- 
ments substantially smaller than the TB model of GSP 
fusf . The LDA, further, leads to a resonant-bond Jahn- 
Teller distortion (as opposed to the usual pairing config- 
uration) for the simple divacancy that is not observed in 
TB calculations. Also, the relaxation energies obtained 
from the GSP-TB model for the divacancies are signifi- 
cantly larger than the corresponding LDA values. 

The formation volumes of the vacancy and the inter- 



stitial have been calculated by Tang et al. |2(J using the 
TB model of Kwon et al. |L2j. The formation volume is 
defined as Af2 = V IC \ ± f2, where V IC \ is the relaxation 
volume associated with the defect (i.e., arising from the 
relaxation of the atoms in the neighbourhood of the de- 
fect) and f2 is the volume per atom of the perfect crystal; 
the plus sign is for vacancies while the minus sign ap- 
plies to interstitials. Using a 216-atom supercell, and 
after a careful search for the equilibrium volume of the 
perfect crystal, Tang et al. obtained a relative formation 
volume of 3% (contraction) for the vacancy and 

— 10% (expansion) for the interstitial. Thus, the volume 
changes arising from the presence of these defects should 
cancel each other to a large extent, in agreement with 
diffuse x-ray scattering experiments [ pll . 

Though the picture is far from being complete, and 
it is therefore difficult to draw meaningful conclusions 
on the accuracy of the TB models, it seems to be the 
case that the model of Kwon et al. 0| overestimates 
the relaxation energies, while that of Lenosky et al. plj] 
appears to be doing better (at least for the monova- 
cancy). Bernstein and Kaxiras p3] have observed that 
the agreement between TB and LDA defect formation en- 
ergies can be improved significantly by relaxing the con- 
straint on the band gap, which is then allowed to vary 
in the fitting process. Clearly, more precise TB models 
are necessary in order to capture the subtle details of 
such low-symmetry situations. Also, more (and better- 
converged) first-principles calculations are needed to pro- 
vide a proper reference database for comparison. 

Rasband et al. J2i| have studied the convergence of in- 
trinsic defect formation energies with respect to potential 
cutoff distance as well as number of points used to sam- 
ple the Brillouin zone. Considering isolated vacancies as 
a test case, they found these variables to affect only very 
slightly the unrelaxed formation energy, while the effect 
of relaxation can be sizable. For instance, increasing the 
cutoff from 3.2 A (between first and second neighbours) 
to 4.1 A (between third and fourth) and using 40 k points 
rather than one causes the vacancy formation energy to 
decrease from 3.67 (cf. Table to 3.15 eV. The corre- 
sponding values for the — , + and 2+ charge states of 
the vacancy are 2.9, 3.6 and 4.1 eV, respectively. The 
+ vacancy is nowhere in the gap a favourable state of 
the defect, and thus gives rise to the so-called "negative- 
W effect, that is an effective correlation energy between 
electrons which is negative [S3] (see also below for GaAs) . 

For the tetrahedral interstitial in its neutral state, now, 
Rasband et al. |2^,^5| find a formation energy of 4.7 eV, 
using a 4.1 A cutoff and 40 k points. This is compara- 
ble to the 4.40 eV value given in Table | (3.2 A cutoff, 
L-point only). In the — , +, 2+, and 3+ charge states, 
the corresponding numbers are 5.5, 4.2, 3.5, and 4.1 eV, 
respectively, taking the Fermi level in the middle of the 
gap. Thus, 2+ interstitials should occur with a much 
larger probability than other charge states. This predic- 
tion of the stability of the 2+ interstitial is in agreement 
with earlier ab initio results and might explain the dis- 
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crepancy between experiment (or more precisely "model- 
fitted" experimental data — cf. Fig. 2 in Ref. |^| ) and 
many calculations. In particular, this result is consistent 
with metal in-diffusion experiments (see p5f for refer- 
ences). 

Rasband et al. [^5| have used TBMD to search for new 
defect structures (interstitials) and have found a whole 
family of them. For the neutral single interstitial, three 
stable configurations are found, viz. the T interstitial 
(formation energy of 4.7 eV), the 110-split interstitial 
(5.0 eV) and the 100-split interstitial (5.4 eV). The cor- 
responding ionised defects with charges between +3 and 
—2 range in energy from 5.0 to 6.3 eV for the 110 split, 
and from 4.5 to 6.3 for the 100 split. As we have seen 
above, the doubly-ionised T interstitial has a formation 
energy of 3.5 eV and is thus much more stable than any of 
the split interstitials. Di-interstitials were also examined. 
Rasband et al. p5| found a novel low-energy configura- 
tion — the "split triple" interstitial, consisting of three Si 
atoms sharing one lattice site and forming an equilateral 
triangle in a (111) plane. The formation energy of this 
defect is a small 3.65 eV per atom in the neutral state, 
dropping to 3.0 eV in the 2+ state. A similar 110 split- 
triple interstitial has a 3.3 eV formation energy, while 
a "Z" configuration has 3.4 eV, both in their 2+ state, 
which is the most stable. There is, to our knowledge, 
no evidence that these objects have been observed ex- 
perimentally, nor are there other calculations to compare 
with. In view of the approximate character of TB, the 
stability of such defects should probably be considered as 
somewhat speculative at this point. 

2. Energy levels 

One advantage of TB over empirical approaches is that 
it gives access to the electronic structure of the material. 
The GSP model does not provide a very good description 
of the band structure of Si. It is nevertheless informative 
to examine, within this model, the effect on the electronic 
properties of the presence of defects. Fig. [I] shows the 
density of states near the band gap for the various defect 
types considered by Song et al. fll7|j . Defects, evidently 
give rise to localised bands in the gap. Table |l| lists the 
positions of the levels associated with the defects in their 
relaxed configurations. 

The LDA results quoted in Table || are those of Ref. 
p6| , where other references can also be found; however, 
only very few calculations of the defect levels have been 
reported. For the monovacancy, for instance, the LDA 
value for the highest occupied level, 0.23 eV, is in rough 
agreement with the self-consistcnt-field calculations of Li- 
pari et al. p6[ , which gives a singlet level at 0.3 eV (but 
relaxation was not taken fully into account). Both cal- 
culations however disagree with the TB result. For the 
divacancies, the error bar on the LDA calculation (about 
0.1 eV) is such that it is not clear that this defect leads 



to levels in the gap, in apparent disagreement with the 
TB results. It should be said, as discussed in Ref. Eg) , 
that the precise positions of defect-induced gap levels 
are rather sensitive to the details of the computational 
model. All calculations, however, point to the fact that 
defect levels move towards the valence-band maximum 
upon relaxing the structure. 

Electron paramagnetic resonance (EPR) experiments 
gives information about the energy levels of charged de- 
fects. The + and 2+ charged states of the vacancy lie 
close to the top of the valence band, at 0.05 and 0.13 eV, 
respectively, whereas the — and 2— lie at 0.57 eV and 
0.11 eV below the bottom of the conduction band. The 
TBMD results of Rasband et al. are in error with experi- 
ment by -0.04, -0.07, -0.17, and -0.90 eV, respectively 
(see [E3[ for experimental references). Except for the 2— 
defect, the agreement can be considered as good. 

3. Vacancy clusters in silicon 

Structural evolution during non-equilibrium processes 
such as irradiation and growth is mediated, to a large 
extent, by the presence of defects or clusters of them. 
Vacancies, for instance, may coalesce and give rise to 
microdefects such as voids, bubbles, or dislocations. It 
is therefore important to have some idea of the structure 
and energetics of such imperfections. 

The case of vacancy clusters has been examined by 
Bongiorno, Colombo, and Diaz de la Rubia (27j] (BCDR) 
using TBMD. This work was motivated to a large ex- 
tent by serious disagreements between model potential 
(Stillinger- Weber) and first-principles studies with re- 
gards, in particular, to the shape and energetics of sta- 
ble clusters. In order to perform these calculations, 
BCDR exploited the power of the Goedecker-Colombo 
linear scaling scheme [ p8| as implemented on a (paral- 
lel) Cray T3E. This allowed very large simulation cells 
to be considered — 1000 atoms — containing vacancy 
clusters varying in size between 1 and 35. In all cases, 
the lowest-energy configurations was obtained through a 
careful relaxation. 

BCDR examined two different cluster 
shapes: spheroidal clusters (SPC), where vacancies are 
created by simply removing atoms from successive radial 
shells about a single vacancy, and "hexagonal ring clus- 
ters" (HRC), where Si atoms are removed following a ring 
pattern in the crystal. These structures evidently are dif- 
ferent from the geometric viewpoint, but also from the 
"chemical" viewpoint. Indeed, for a given cluster size, 
SPC structures have more dangling bonds than HRC 
structures, as demonstrated in the top panel of Fig. |^, 
and SPC are therefore expected to be less favourable than 
HRC because dangling bonds are expensive. This is in 
fact confirmed by the formation energy, middle panel in 
Fig. ||, at least for clusters containing 1-24 atoms. For 
larger clusters, the aggregation path is likely very com- 
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plex and probably depends on the details of the kinetics 
of the non-equilibrium process (irradiation or growth). 

The binding energy for the various clusters is given 
in the bottom panel of Fig. ||; there are clearly some 
"magic clusters", i.e., clusters who are much more stable 
than others slightly different in size. This is the case, for 
instance, of HRC clusters with n = 6, 8, 10, 12, 14, 16, 
18, etc.; these magic clusters are the result of minimising 
the number of dangling bonds as as well as structural 
rearrangements (internal reconstructions). 



B. Impurity levels in Si and GaP 

As we have seen, relaxation affects strongly the struc- 
ture and energetics of intrinsic defects in semiconductors. 
This is true also of impurities. The deep levels associated 
with several impurities in Si and GaP have been investi- 
gated using TBMD by Li and Myles f||]. Tneir approach 
is based on the theory of Hjalmarson et al. |3(| of deep 
levels, augmented to include lattice relaxation by adding 
a repulsive pair potential derived from Harrison's over- 
lap interactions |31j . The host material is described using 
the sp 3 s* TB model of Vogl et al. jj2j. The relaxation is 
performed via MD, but restricted to nearest neighbours 
and Td symmetry-conserving displacements, i.e., breath- 
ing modes. Full details of the method can be found in 
Ref. M. 



Table III gives the energies of several deep-level impu- 
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rities in Si and GaP as calculated by Li and Myles 
These are of A\ (as well as Td) symmetry, i.e., s-like. 
Also given in Table [II are values from experiment (see 
p9[ for references) and the results from the Hjalmarson 
et aZ.'s theory [pQ] , which does not include relaxation. 
Clearly, relaxation is sizable in all cases and improves 
quite significantly the agreement with experiment. Re- 
laxation proceeds inwards in all cases except GaP:0. Li 
and Myles indicate that inclusion of second-neighbour 
relaxation changes the results very little. 



C. Boron in silicon 

Boron is a dopant which is routinely used in the semi- 
conductor fabrication process and it is therefore impor- 
tant to understand at the fundamental level how the host 
material is affected by the impurity, and how the lat- 
ter diffuses. As a first step in this direction, Rasband 
et al. (HJH] have used TBMD to study defect-dopant 
pairs as well as boron interstitials in silicon. For Si-Si 
interactions, the GSP model was employed. For Si-B in- 
teractions, a new model a la GSP was developed, with 
the parameters determined by fitting to ab initio band- 
structure energies for zinc-blende SiB; full details can be 
found in Ref. g|. 

The TBMD results are found to be generally in good 
agreement with the ab initio calculations of Nichols et al. 



p3| . For the boron T interstitial, the TB model gives a 
formation energy of 3.7 eV, compared to 3.9 eV from first 
principles. For the neutral boron-substitutional-vacancy 
complex, Rasband et al. find 2.9 eV compared to 3.0 for 
Nichols et al.; the binding energy of this complex is —0.39 
eV, which compares well to —0.5 ab initio. The bind- 
ing energy of the boron-substitutional-interstitial com- 
plex is —1.2 eV from TBMD, in line with model-fitted 
data, —1.5 eV. Relaxation of the four Si atoms near a 
B substitutional impurity is found to be small — about 
0.04 A, independent of the state of charge. 

Rasband et al. |25| have looked for "new" defects in- 
volving boron. Starting with a boron atom in a T- 
interstitial position with charge states in the range —2 
to +2, they identified a series of seven different inter- 
stitial species having formation energies, in their most 
stable charge states, between 3.1 and 5.5 eV. They also 
examined the case where the starting structure was a 
combination of Si 110-split interstitial adjacent to a B 
T-interstitial; the resulting configuration after relaxation 
was a 110-split di-interstitial with a formation energy of 
2.5 eV per atom. Further, they observed the formation 
of a "ring" di-interstitial having a formation energy of 
3.1 eV per atom exhibiting a distinctive five-membered 
ring. 

The above results, just as was the case for native 
defects in silicon, indicate that new elements must be 
taken into account when experimental defect concentra- 
tion data are analysed. In particular, assumptions re- 
garding the prevailing charge state of any particular de- 
fect should probably be re-examined at the light of the 
above (and other) results. 



D. Point defects in GaAs 

In the above study of impurity levels in Si and GaP, 
Section B, relaxation was restricted to radial displace- 
ments. While such distortions might constitute the most 
important component of relaxation, it may be expected 
that, quite generally, the Td symmetry is broken by the 
presence of defects, and this can only be revealed through 
full relaxation of the lattice. 

A TBMD investigation of defect relaxation and ener- 
getics in GaAs has been proposed by Seong and Lewis 
p4j . The study is based on the model of Molteni, 
Colombo, and Miglio [jj5|, which takes charge-transfer 
effects into account through a screened, and thus short- 
range, Coulomb term. Native point defects in GaAs play 
an important role; this is true, in particular, of the so- 
called "EL2" defect, which can compensate residual ac- 
ceptor impurities and pin the Fermi level at midgap, so 
that GaAs can be manufactured without intentional im- 
purity doping. This defect is thought to be related to the 
As antisite (substitution of a Ga by an As). 

For tetragonally-distorted configurations, it is conve- 
nient to describe the relaxation in terms of three dis- 
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placement components [p[|3q| : the "breathing" (radial) 
mode, and two "pairing" modes, which measure the lat- 
eral displacements of the ions (i.e., perpendicular to the 
breathing mode). Table |^ gives the breathing- mode dis- 
placements for Ga and As vacancies and antisites in var- 
ious states of charge. The pairing modes are negligibly 
small for Ga vacancies and As antisites, but are sizable 
for As vacancies and Ga antisites; they are given in Table 

Ga vacancies are found to relax inwards by an amount 
which is independent of the state of charge, leading to an 
open- volume contraction of about 34%, consistent with 
positron lifetime measurements |37|,[58j, as well as first- 
principles calculations p6| , though relaxation in the lat- 
ter case is substantially smaller and the position of the 
corresponding localised states different. The relaxation 
of As antisites also is the same for all charge states; it is 
small and, in contrast to Ga vacancies, proceeds out- 
wards, leading to a volume expansion of about 7.5%. 
This is also what is found in first-principles calculations 
p9| fm an d experiment We find the localised state 
for this defect to lie 0.49 eV above the valence band max- 
imum, in accord with tunneling spectroscopy measure- 
ments @. 

Both As vacancies and Ga antisites undergo 
tctrahcdral-symmctry-breaking distortions, i.e., have 
non-zero pairing displacements. In both cases, further, 
the relaxation pattern depends sensitively on the state 
of charge, as can be seen in Tables |fv] and [v|. For As 
vacancies, changes in the open volume range from a 54% 
expansion in the + state to a small 3% contraction in the 
2— state. The positron lifetime is therefore expected to 
be longest in the + state and shortest in the 2— state; 
this is consistent with experiment, which finds a decrease 
from 295 to 258 ps when the charge changes from neu- 
tral to singly negative p7| . In the case of Ga antisites, 
the changes in the open volume vary between +50% and 
a negligible —1%; the latter is for the doubly charged 
vacancy and corresponds, approximately, to the GaAs 
crystal in its "normal" state. (An As with five valence 
c~ is replaced by a Ga with three e _ , then two electrons 
are added to it). 

It is of interest to examine more closely the structure 
of the Ga antisite. The singly-negative Ga antisite un- 
dergoes a small (<2%) outward breathing-mode displace- 
ment, and similar pairing-mode relaxations of the nearest 
neighbours. The antisite itself is displaced slightly from 
its ideal lattice position [by about 0.1 A in the (111) di- 
rection]. The distance between the antisite and one of 
its Ga neighbours (referred to as "Gai" by Zhang and 
Chadi ^3|) is 2.65 A, a little bit longer than the equilib- 
rium bond length (2.45 A). The distance between the Ga 
antisite and the three other Ga neighbours, which relax 
very little, is 2.44 A. This atomic structure near Ga^ s is 
very similar to that of Zhang and Chadi |43) who studied 
the problem from first principles; the configuration is de- 
picted in Fig. ||(a). The relaxed state of the neutral Ga 
antisite is shown in Fig. ||(b). The defect moves away by 



0.4 A from its equilibrium position; the Gai atom also 
moves away, but in the opposite direction, by 0.48 A. This 
leads to a "broken-bond" configuration with a relaxed 
Ga^ s -Gai distance of 3.32 A, in perfect agreement with 
Zhang and Chadi's calculation Ea]. Positively-charged 
Ga antisites show similar relaxations. 

The variation with electron chemical potential (fi e , 
where < /i e < E g and E g is the width of the gap) of 
the formation energy of a given defect in various states 
of charge provides information about its ionisation level 
and most stable state. This is displayed in Fig. |] for the 
four types of defects considered. A detailed discussion 
of these results can be found in [[34J. The formation en- 
ergy also depends, however, on the atomic environment 
that prevails during growth, i.e., on atomic chemical po- 
tentials. Fig. U shows, for three different values of /i e , 
the defects of lowest energy as a function of the devia- 
tion A/iQa from the bulk value. At both the valence-band 
maximum (VBM) and midgap, the As antisite is the most 
favourable defect in the As-rich limit, while the Ga an- 
tisite is the one preferred in the Ga-rich limit. In the 
middle of the A^iq^ range, the vacancies V£ s and Vq~ 
take over at midgap, though, in only a narrow window 
of values. At the conduction-band maximum (CBM), 
now, Vq~ dominates in As-rich conditions, while Ga^~ 
is favoured in the Ga-rich limit. Thus, overall and by in 
large, antisites have the lowest energies, and these defects 
are expected to be much more prevalent than vacancies 
in GaAs. This is consistent with the current interpre- 
tation fHHhH] of the EL2 defect being a neutral As 
antisite with tetrahedral symmetry. 



E. Dynamics and kinetics of point defects 

1. Vacancy and interstitial diffusion and recombination 

Vacancies and interstials (as well as small complexes 
of them and impurities) contribute strongly to the trans- 
port of mass in crystalline solids. Yet, little is known on 
their equilibrium concentrations and diffusivities. The 
values from experiment vary widely, and it is not clear 
which of vacancies or interstials dominate transport at a 
given temperature. Further, the recombination process 
of the two species is not well understood, and in partic- 
ular it is not clear whether or not there exists an energy 
barrier against the recombination. These questions have 
been addressed, in part, using MP and empirical poten- 
tials (e.g., Stillinger- Weber jl{|); such calculations how- 
ever cannot account for the subtle quantum-mechanical 
details of the process. First-principles methods, on the 
other hand, are hampered by limitations in size and time 
that make the study of such problems almost impossible 
at this time. TBMD offers a good compromise between 
the two, and was used by Tang et al. to study vacancy 
and interstitial diffusion and recombination |20|] . For this 
problem, they employed the TBMD model of Kwon et al. 
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p2| , which appears to be more accurate than the GSP 
model insofar as the properties of point defects are con- 
cerned for a 216-atom model. The diffusivity calculations 
were however performed on a 64-atom (plus or minus one) 
system, running for as long as 200 ps. 

The diffusion constant is given by the product of diffu- 
sivity and equilibrium concentration, D(T) = d(T)C(T), 
where 



and 



d{T) = d Q cxp{-E m /k B T) 



C(T)=e X p[-(Ef-TS f )/k B T]. 



E m is the migration energy and Ef is the formation en- 
ergy. If the entropy of formation, Sf, is assumed to be 
independent of temperature, then the diffusion constant 
takes the Arrhenius form: 

D(T) = A, exp{-E a /k B T) 

where Dq = d Q exp(Sf /k B ) and E a = E m + Ef. D(T) is 
indeed observed in crystalline silicon to be Arrhenius over 
a wide range of temperatures. The precise values of the 
prefactor Dq and activation energy E a for the individual 
species are still a matter of discussion, but the follow- 
ing values have recently been given by Gosele et al. [pfj : 
E% = 4.03 eV (vacancies), E* a = 4.84 eV (interstitials), 
and BljDl = 1.53 x 10 4 . 

The diffusivities of vacancies and interstitials were sim- 
ulated explicitly by Tang et al. Q using TBMD. The 
migration energies were found to be E^ = 0.1 eV and 
.E^j = 1-37 eV. In this model, the formation energies are 
Ej = 3.97 eV and Ej = 3.80 eV, leading to activation 
barriers of E^ = 4.07 and E^ = 5.17 eV, in remarkable 
agreement with the experimental values quoted above. 
Tang et al. did not calculate the entropies of formation 
- this is a very difficult calculation — but have fitted 
them to first-principles calculations and experiment; they 
obtain Sj = 9k B and Sj = 11.2k B . The resulting pref- 
actors are Dq = 0.96 cm 2 /s and = 1.16 x 10 4 cm 2 /s, 
i.e., Dq/D^ ~ 1.21 x 10 4 , again in agreement with ex- 
periment. Arrhenius plots of the diffusion constants are 
given in Fig. ^|. Diffusion is dominated by vacancies at 
low temperature. At high temperature, above 1080 °C, 
in spite of a sizeably larger migration barrier, interstitials 
take over. This crossover effect is due to the relative val- 
ues of the prefactors — Dq is four orders of magnitude 
larger than Dq — and find it's origin in the Meyer-Neldel 
(compensation) law, which states that, for a family of ac- 
tivated processes, the prefactor increases exponentially 
with the activation barrier. The validity (and origin) of 
the Meyer-Neldel law was established in the case of sur- 
face diffusion by Boisvert et al. . 

Mass transport, as noted above, will evidently be af- 
fected by the rate at which vacancies and interstitials re- 
combine, which itself depends on the existence or not of a 
recombination barrier. Tang et al. 20 have used TBMD 



to investigate this problem. Starting with an interstitial- 
vacancy pair with the two defects separated by some dis- 
tance — between two and six nearest-neighbour spacings 
- MD simulations were carried out at finite tempera- 
ture in order to follow the migration path and eventual 
recombination of the defects. The latter does not al- 
ways take place. For instance, using as initial configura- 
tion a < 110 > dumbbell interstitial and a vacancy three 
nearest-neighbour spacings away, Tang et al. found that 
recombination did not occur at 300 K: thermal activ- 
ity at this temperature is not large enough to overcome 
the local distortions around the dumbbell induced by the 
close approach of the vacancy, at least on the timescale 
of the simulations. Instead, an "I-V complex" forms; the 
quenched-in state of the complex, which has a formation 
energy of 3.51 eV, is shown in Fig. 0(a). At elevated tem- 
perature, however, the I-V complex annihilates. This is 
shown in Fig. ^(b)-(e). The annihilation path is found to 
be a bond-switching process and takes place in about 10 
ps at 1500 K; the activation (migration) barrier E^ for the 
process is about 1.23 eV. The lifetime r of the complex is 
roughly given by r _1 = v exp(— Eb/k B T); taking v as the 
Debye frequency, ~ 10 13 Hz, one finds r to be of the order 
of hours at room temperature, but only a few microsec- 
onds at typical annealing temperatures (i.e., 700-800 K 
or so). The recombination of vacancy-interstitial pairs, 
in particular with regards to the electronic structure, was 
examined by Cargnoni et al. [^2 53 1 by a combination of 
TBMD and ah initio Hartree-Fock calculations. 



2. Hydrogen diffusion in silicon 

Because it may form complexes with a variety of in- 
trinsic or extrinsic defects, hydrogen in semiconductors 
affect deeply their optical and electronic properties. It 
is usually present as a result of the complicated fabri- 
cation process, but is also often intentionally included 
in order to passivate defects (e.g., unsaturated, or "dan- 
gling", bonds). Being a light species, hydrogen diffuses 
readily, inducing additional defects along the way and 
thus affecting the transport properties of the material to 
an extent which is determined by its relative concentra- 
tion. It is therefore important to understand diffusion at 
the atomic level in order to gain better control on the 
properties of semiconductors. 

The rate of diffusion of hydrogen in crystalline sili- 
con remains a matter of debate in spite of the simplicity 
of the structure of the host material. At low tempera- 
tures (below 800 K), the experimental values that have 
been reported in the literature vary widely, sometimes 
by as much as two orders of magnitude at a given tem- 
perature (see, e.g., 54 for references) . There are very 
few data available at high temperatures (above 1000 K), 
but they are in reasonable agreement with the ah initio 
MD simulations of Buda et al. |35) for H + (proton) diffu- 
sion. Unfortunately, because the timescale for diffusion 
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increases exponentially with temperature, ab initio sim- 
ulations cannot be carried out at lower temperatures and 
there exists no empirical model that can deal with this 
system in a sufficiently accurate way. 

The problem of hydrogen diffusion in silicon was ad- 
dressed using TBMD by Panzarini and Colombo (PC) 
|| as well as by Boucher and DeLeo (BDL) @. This 
constitutes an interesting application of the method as it 
extends quite significantly the range of temperatures that 
were covered by the ab initio simulations of Buda et al. 
p5| while accounting, still, for the quantum-mechanical 
nature of the system. For Si-Si interactions, both used 
the GSP model; for other interactions, the appropriate 
parameters were determined by fitting to selected prop- 
erties of the silane (SiEL) molecule. 

Both models give the bond-center (BC) site as the 
equilibrium site in static conditions (i.e., at K), in 
agreement with electron paramagnetic resonance and 
muon spin resonance experiments |5S|] , as well as previ- 
ous first-principles calculations (see, e.g., |||]). In the 
BDL model, this accord is obtained at the expense of 
a phenomenological parameter that accounts for differ- 
ences in environment between silane and crystalline Si. 
Other equilibrium sites, almost degenerate in energy (dif- 
ferences of less than ~0.1 eV) with the BC site, are also 
found. PC, for instance, find the hexagonal interstitial 
site to be very close in energy to the BC site. Both cal- 
culations predict that when H sits in the BC position, 
nearby Si atoms relax by about 0.4 A, while very little 
relaxation is observed when H is in the hexagonal site; 
this agrees with first-principles calculations. 

The diffusion path is also subject to debate. Though 
experiment and theory agree that the lowest energy site 
is the BC site, the energy of nearby metastable sites — 
which likely play an important role in diffusion — is not 
known precisely. Details of the diffusion path, further, 
are expected to be strongly affected by dynamical effects 
because the heavy Si atoms cannot follow adiabatically 
the motion of the lighter H atom. The ab initio simula- 



tions of Buda et al. 55 1 , which cover a timescale of about 



4 ps, suggest that diffusion consists of a series of hops be- 
tween highly-symmetric interstitial sites, but other pos- 
sibilities can also be envisaged. 

Both PC and BDL carried out their TBMD simulations 
for a single H in a 64-atom c-Si supercell. The simulations 
by BDL cover the temperature range 1050-2000 K, and 
run for a maximum of 42 ps, while PC examined temper- 
atures in the range 800-1800 K, running their simulations 
for as long as 300 ps at the lowest temperatures. 

Fig. |^ presents the results of several measurements 
of the diffusion constant, plotted in the manner of Ar- 
rhenius. Also indicated on this plot are the ab initio 
MD data points of Buda et al. pEfl ; they are found to 
agree (at least qualitatively) with the high-temperature 
experimental points of Van Wieringen and Warmholtz, 
which can be fitted to an Arrhenius law, D(T) — 
D exp(-E A /k B T), with D = 9.41 x 10~ 3 cm 2 /s and 
Ea = 0.48 eV. When extended to low temperatures 



(dashed line in Fig. g), one clearly sees the deviations 
from the Arrhenius behaviour; it should be said, how- 
ever, that there is no "guarantee" that diffusion should 
be Arrhenius over the whole range of temperatures. 

The TBMD data of PC are also indicated in Fig. |. 
The agreement with experiment is clearly excellent in the 
high temperature limit. The data of BDL are not plot- 
ted in this figure, but they are found to be extremely well 
fitted by the Arrhenius law with Dq = 6.91 x 10~ 3 cm 2 /s 
and Ea = 0.45 eV all the way down to 1050 K. This is 
in striking agreement with experiment, as can be judged 
by the close similarity between the prefactors and energy 
barriers. (The differences are unsignificant). PC how- 
ever observe deviations from Arrhenius already at 1200 
K. The reason for this minor discrepancy with the BDL 
data can perhaps be traced down to statistics. Indeed, 
300 ps remains relatively short on diffusion timescales. 
For instance, for a diffusion constant of 10~ 6 cm 2 /s, a 
quick calculation indicates that the average distance vis- 
ited by a diffusing particle would be about 4 A. This 
might explain, in part, the deviations that are seen at 
even lower temperatures; clearly, however, the TBMD 
data are consistent with the low-temperature behaviour 
seen in experiment. Further calculations are evidently 
necessary to reconcile theory and experiment at low tem- 
peratures. 

Detailed analysis of the trajectories of individual atoms 
makes it possible to elucidate the mechanism for diffu- 
sion. Fig. ^ shows the diffusion path of the H atom over 
a 25-ps period at 1200 K as calculated by PC; BDL obtain 
very similar results. The H atom diffuses preferentially 
via a sequence of jumps from one BC site to another, 
spending very little time in between; this corresponds to a 
low-memory, high-friction regime, where the directions of 
consecutive jumps are uncorrelated. It is found, further, 
that the H atoms avoids the low charge-density regions, 
which is inconsistent with the observations of Buda et 
al. H for the diffusion of H+. Both PC and BDL find 
a correlation between diffusion events (jumps) and the 
vibrations of nearby Si atoms: At low temperature (less 
than 850 K or so), the Si-Si stretching mode (~65 meV) 
is not thermally excited, and diffusion becomes difficult 
because the BC site is an equilibrium site for hydrogen. 
At these temperatures, further, long jumps (i.e., to sites 
more distant than near-neighbours) are quenched in. 



IV. DISORDERED PHASES 

As already noted in the Introduction, disordered ma- 
terials — covalently-bonded semiconductors and chalco- 
genides — is another area where TBMD simulations have 
provided important new insights, in particular with re- 
gards to the interplay between structure and electronic 
properties. 

Many empirical potentials have been developed that 
give a reasonable description of the liquid and crystalline 
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phases of semiconductors, as well as, for some, clusters; 
this is the case for instance of the models of Stillinger 
and Weber ]49| , Tersoff ]5{|, Biswas and Hamann ]6C| ] 
and, recently proposed, Bazant et al. |6l|] These mod- 
els however lack the generality to account for the large 
variety of highly strained environments — caused by elas- 
tic, topological or chemical disorder — that are encoun- 
tered in amorphous materials. Thus, for instance, the 
Stillinger- Weber potential fails to reproduce the clean 
separation observed experimentally between first and sec- 
ond neighbour peaks in the radial distribution function. 
The problem is particularly serious for materials such as 
amorphous carbon, where atoms can be found in several 
different bonding states, or compound materials, where 
bonding is partly ionic and the nature of the chemical 
environment plays an important role. There exists, for 
example, no satisfactory empirical potential for the III-V 
amorphous compounds. 

The availability of accurate potentials is a necessary, 
but not sufficient, condition for constructing structurally 
sound amorphous models. The "preparation" algorithm 
must be capable of finding a reasonable ground state — 
one in which the number of topological and electronic 
defects is a minimum — i.e., relaxation must be as ex- 
haustive as possible. Unfortunately, this is more difficult 
with TB than empirical potentials. 

The strong intermediate range disorder found in amor- 
phous semiconductors poses an additional challenge for 
the development of semi-empirical TB potentials. As 
for the liquid phase, near-neighbour shells are often ill- 
defined: The first and second neighbour shells can be 
blurred either through disorder or simply large mismatch, 
as is the case for InP, for example. Attention must there- 
fore be given to the radial behaviour of the empirical 
repulsive potential as well as the cutoff of the matrix el- 
ement interactions. In spite of these concerns, current 
semi-empirical potentials, even though often fairly short- 
range, yield good agreement with experiment. Moreover, 
ab initio TB interactions, which normally include longer- 
range interactions, are not affected by this situation. 

The first amorphous materials to be studied with 
TBMD were the elemental semiconductors silicon and 
carbon. Silicon, which has a well-defined sp3 bonding, is 
a relatively straightforward choice; TBMD a-Si models 
are discussed in the next section. Carbon, on the other 
hand, is much more difficult because of the added diffi- 
culty associated with its many bonding states; another 
article in this special issue is dedicated to a-C and we 
therefore do not discuss it in any detail here but just 
give, for completeness, a brief overview of work on this 
material. 

Hydrogenation of amorphous semiconductors reduces 
significantly the strain and greatly improves the elec- 
tronic properties by removing deep states in the gap. 
This is particularly so for a-Si, which requires consid- 
erable amounts of H to achieve device-quality electronic 
properties. TBMD simulations have been run to study 
its effect on structural and electronic properties, in spite 



of problems associated with its small mass and high dif- 
fusivity. Likewise, TBMD has been used to investigate 
compounds (mostly binaries) — GeSe2, GaN, GaAs - 
where additional complications arise from partly ionic in- 
teractions. 



A. Elemental semiconductors 

1. Amorphous silicon 

Silicon is the material of reference in the study of amor- 
phous semiconductors and has been simulated numeri- 
cally using a variety of techniques. It is widely accepted 
as a realisation of Polk's idealised continuous random 
network (CRN) (62|, which regards the material as a 
collage of randomly-oriented, corner-sharing tetrahedra, 
thus possessing perfect coordination. Algorithms have 
been devised for constructing Polk-type CRN's on the 
computer, and these will serve as a reference for simula- 
tions based on total-energy minimisation such as TBMD. 

As noted earlier, amorphous samples can be prepared 
in the computer in many different ways. Quenching from 
the melt and annealing is a popular method because it 
is akin to the real fabrication process. Quench rates 
for computer models (in particular TBMD) are how- 
ever many orders of magnitude larger than real ones 
and the main difficulty therefore resides in cooling slowly 
enough for the system to be able to find a reasonable 
low-energy state. Monte-Carlo methods such as the 
Wooten- Winer- Weaire bond-exchange algorithm |33 64 
and the activation-relaxation technique of Barkema and 
Mousseau [ ]65| , |S6[ , can yield lower-energy configurations, 
but at the expense of a heavier computational effort un- 
less they are conducted, in part, using empirical poten- 
tials. 

Kim and Lee @ (KL) have used TBMD (with the 
GSP model) and the melt-and-quench approach to con- 
struct a 64-atom model for a-Si. The liquid was produced 
by running at high temperature (1750 K) for 8 ps. (The 
liquid phase itself is much better described by the GSP 
TB interaction than by the SW potential; in particular, 
the number of nearest neighbours and the angular distri- 
butions are in close agreement with the values obtained 
from ab initio MD |68|j70[| ) . The cell was then cooled 
at a rate which amounts to 10 15 K/s. The simulations 
were done in the microcanonical ensemble at the density 
of the liquid, about 10% larger than that of the crystal, 
which is itself 1.6% denser than the amorphous phase 
jnj. Thus, the final amorphous phase is under high com- 
pressive strain. The structural properties of the result- 



ing structure are given in Table VI, where they are also 



compared with the predictions of other models discussed 
below. The total coordination number, 4.28, is large, a 
consequence of the excessive cooling rate but also of the 
compressive strain on the system which severly hampers 
relaxation in view of the short period of time covered by 
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the simulation. The model displays a rather clean sepa- 
ration between first- and second-neighbour shells despite 
an unrealistically large number of defect states (float- 
ing bonds) in the electronic bandgap, mostly due to the 
strain associated with overcoordination. 

Servalli and Colombo [[72| have studied the influence of 
the cooling rate on the structure the material. They con- 
sidered a 64-atom system, first ran for 27 ps in the liquid 
phase, then cooled at rates in the range 0.096 — 4 x 10 14 
K/s — that is up to 100 times slower than in KL's sim- 
ulations. The volume was varied linearly with tempera- 
ture between values appropriate for the liquid and amor- 
phous phases. Overall, the coordination is in much bet- 
ter agreement with experiment than KL's model. More- 
over, there appears to be a strong correlation between the 
cooling rate and the "quality" of the sample. From vi- 
sual inspection of the data plotted in the original article, 
the bond angle distribution, in particular, decreases by 
about 20% between fastest and slowest cooling rate, com- 
ing quite close in the latter case to experimental values. 
Even the slowest rate, however, might not be sufficient 
for proper relaxation: In one run, for instance, a "major" 
relaxation event was observed after about 57 ps, indicat- 
ing that the relaxation time is much longer than a few 
phononic periods even in such small systems. Servalli and 
Colombo also prepared a 216-atom a-Si model [[72] which, 
because of computational limitations, was quenched-in 
over a fairly short time of 11.4 ps. In spite of this, the 
electronic density of states is comparable to that obtained 
from a smaller system relaxed for a longer time, suggest- 
ing that the size of the unit cell might be an important 
factor in the total stress found in computer-generated 
samples. 

In order to minimise the computational workload im- 
posed by the long quenching process, mixed approaches, 
using empirical potentials for the time-consuming prepa- 
ration followed by TBMD relaxation, can also be used. 
Mousseau and Lewis |7^,|7^], for example, have combined 
TBMD with an efficient optimisation scheme based on 
classical potentials. A "rough" 216-atom model was first 
prepared from a randomly-packed configuration using the 
activation-relaxation technique ]65| , |66[ | together with a 
Stillinger- Weber-type potential. The resulting structure 
was then relaxed using TBMD and the GSP potential. 
The final TBMD stage ensures that the model is physi- 
cally realistic; the initial empirical-potential relaxation is 
therefore, to some extent, artificial — as is also the case 
of the Wooten-Weiner-Weaire bond-switching process — 
but ensures, when combined with such a powerful opti- 
misation scheme as the activation-relaxation technique, 
that the structural model is fully optimised. This is actu- 
ally demonstrated in Table VI : the average coordination 
number is very close to four (but slightly below), as ex- 
pected, and perhaps more significant, the width of the 
bond angle distribution is much smaller than obtained in 
any other model. 

An alternative approach was proposed by Drabold et 
al. W^M; it consists in "incompletely melting" the sample 



before cooling. This is achieved by introducing a vacancy, 
which destabilizes the lattice: because the system is small 
(about 64 atoms in this case), the melting temperature 
decreases markedly and relaxation proceeds more easily. 
After a short microcanonical run at high temperature, 
the system is taken down to zero temperature and re- 
laxed; the calculations were performed using the Sankey- 
Drabold "06 initio TB" scheme |7q| . Although this tech- 
nique allows for a rapid production of samples, it does 
not seem at present to be able to yield quality amor- 
phous networks; the radial distribution functions (RDF) 
of the samples prepared in Ref. |75| either contain traces 
of crystallinity after quench or show the presence of high 
levels of strain. A more controlled procedure could, in 
principle, lead to relatively good structures without re- 
quiring full melting of the crystal. 

The atomic structure of a-Si is not known in detail 
from experiment and it is therefore difficult to assess the 
various computer models. Computer model preparation 
is therefore a challenge by itself as much as a necessary 
first step for further studies. For these reasons, relatively 
little effort has been spent on the actual properties of 
the material. One exception is the study by Drabold, 
Fedders and co-workers of dynamical fluctuations, both 
structural and electronic [Q. Using a 63-atom unit cell 
obtained using the method discussed above, they followed 
the electronic states as a function of time and showed that 
these could fluctuate significantly even at 350 K |f77| . In 
particular, they observed the lowest unoccupied molecu- 
lar orbital to show larger fluctuations than the states at 
the top of the valence band. Moreover, localised states 
fluctuate more than the extended states. Although the 
latter observation can be understood simply in geomet- 
rical terms — the effective mass of the localised state is 
smaller than that of the extended state, the former one 
is not well explained and requires more detailed char- 
acterisation. Studying the structural consequences of 
charged defects in their unit cell, Fedders et al. found 
that an electronic transition can induce structural rear- 
rangements that involve up to many tens of atoms J78| . 
The size of the unit cell prevents quantitative predictions 
to be made, but these results raise interesting questions 
which need to be addressed on larger networks. 

Another application is that of De Sandre et al. [79| ] 
who have computed the elastic constants as a function 
of temperature for the 64-atom TBMD model prepared 
by Servalli and Colombo (72). At finite temperature, the 
elastic constants arc defined as the sum of three contri- 
butions: potential, kinetic and fluctuations Jf9| . Com- 
parison with experiments and results from the empirical 
SW model reveals that the TB potential does better than 
SW but some discrepancy with experiment remains. As 
a general trend, the elastic constants seem to soften as a 
function of temperature. The exact relation, however, is 
somewhat hindered by significant fluctuations that could 
be related to slow relaxation processes taking place dur- 
ing the simulation. 

The physical origin and the density of defect states in 
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the bandgap of a-Si — in particular the band tails — is 
of both fundamental and practical interest. Experimen- 
tally, it is known that a-Si contains up to 1% of defects 
|fnj| , preventing its use in electronic devices. Direct com- 
parison with experiment requires large unit cells in order 
to provide a proper description of the gap region, and 
thus TBMD cannot be used for this purpose. Large em- 
pirical models can however be used "as is" in order to 
compute the TB electronic structure, since this requires 
a single matrix diagonalisation. This was done by Mer- 
cer and Chou as well as Holender and Morgan fjO[ , 
who examined 588- and 13824-atom unit cells, respec- 
tively. Figure 10 shows the electronic density of state for 



a series of models relaxed with either the Keating or the 
Stillinger- Weber potential . The density of states in 
the gap clearly correlates with the number of coordina- 
tion defects. As shown in Fig. |l^ (b) (see also Ref. @), 
there is, however, no one-to-one correlation between co- 
ordination and electronic defects. Some coordination de- 
fects yield states which are deep in the valence band while 
four-fold coordinated atoms with highly-distorted envi- 
ronments can form trap states at midgap. The "rule" 
that arises is therefore that highly-strained networks, 
even perfectly coordinated, are more prone to give rise to 
defects in the electronic gap than low-energy amorphous 
structures with even a few coordination defects. 



Recently, Dong and Drabold (DD) |81| have reported 
a detailed study of the band-tail states in an unrelaxed 
4096-atom Wooten- Winer- Weaire model of a-Si produced 
by Djorjevic et al. |p^| . DD find that the valence band 
tail is well described by an exponential decay p(E) = 
exp(—E/E ) with E = 190 meV. These band-tail and 
gap states also show a significant degree of spatial local- 
isation (Fig. [ll]). This implies that these states cannot 
conduct current unless their density is such that there is a 
significant overlap between them. Showing that many lo- 
calised states actually do overlap significantly with other 
states of similar energy, DD proposed a mechanism - 
the "resonant-cluster proliferation" — that could lead to 
conduction by percolation of overlapping localised states 
of similar eigen energy through the whole system, thus 
explaining the existence of a "mobility edge" in amor- 
phous materials. 



activity in this field. 

The structure of a-C depends on the method of prepa- 
ration: material prepared by evaporation or sputtering 
tends to be sp2 rich while that produced by mass-selected 
ion-beam deposition presents large concentration of sp3 
(diamond-like) carbon. This situation is also found nu- 
merically: Comparing results for 216-atom unit cells of 
a-C prepared at different densities, Wang et al. |83|] found 
that the static structure factor of the low density cell (2.2 
g/cm 3 ) gives best agreement with experimental data for 
sputtered a-C [ p4[ . In this structure, prepared from the 
liquid phase by quenching at a rate of 5 x 10 14 K/s, 80.6% 
of the atoms are three-fold coordinated, while 7.4 and 
12% have four and two neighbours, respectively. Wang 
et al. further observed that the large density of three- 
fold carbons, which form compact regions surrounded by 
two- and four- fold coordinated atoms, leads to fully con- 
ducting — and electronically uninteresting — materials. 
Much effort has therefore been directed towards generat- 
ing dense, and insulating, tctrahcdral amorphous carbon 
using different TB models |8^-|8^|, as we discuss next. 

Lee et al. have developed a semi-empirical TBMD 
model which appears to provide a reasonable description 
of the structural properties of tetrahedral a-C. However, 
the complete absence of a gap in the electronic density of 
states reveals some inadequacies in the model. Also us- 
ing a semi-empirical TBMD model, Wang et al. |^5| have 
studied dense a-C and found results consistent with the 
predictions from the more accurate " ab initio- type TB" 
models of Drabold et al. ^36) and Kohler et al. [p8[ . A 
most surprising feature of these models is the presence of 
a single wide band in the vibrational density of states in 
the range 300-1400 cm -1 that contrasts markedly with 
the sharp structures found in both graphite and diamond. 
As the density of the material decreases, the band splits 
into two wide peaks. Fig. |l2| shows the vibrational den- 
sity of state for 128-atom models of a-C at different den- 
sities |Q. It has been proposed by Kohler et al. |8g| l 
that this peak structure arises from the large strain vari- 
ations about the C tetrahedra, but this remains to be 
established more precisely. 



B. Hydrogenated amorphous silicon 



2. Amorphous carbon 

TB studies of carbon are the object of another article 
in this Special Issue; we therefore do not provide here an 
exhaustive review. For completeness, however, we men- 
tion some relevant work. The modelling of amorphous 
carbon is complicated by the many bonding states that 
carbon can exist in — sp, sp2 and sp3 — and devel- 
oping a TB potential that can properly account for the 
delicate balance between these three states is difficult. 
Nevertheless, because of the fundamental and technolog- 
ical importance of the material, there has been a lot of 



Because it contains a high concentration of defects, 
which depends strongly on the mode of preparation, a- 
Si usually cannot be used dir ectly in electronic devices. 
As already noted in Section HIE 2, hydrogen is often 
intentionally included in the material so as to saturate 
the dangling bonds, thus reducing the concentration of 
defects to an acceptable level. Proper models of a-Si 
are therefore prerequisite to a detailed study of a-Si:H. 
Describing the interactions of H with Si within a TB 
scheme however is a difficult exercise, not only because 
of the added complexity of multi-atom interactions, but 
also because of the subtle quantum-mechanical bonding 
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properties of hydrogen. Several models have been devel- 
oped over the last few years 54 , 5(],^|h| , based on either 
the semi-empirical GSP scheme or the ab initio Sankey- 
Drabold Harris-functional approach (SD-TB) |7q| ; the 
latter defines in a more natural way the long-range inter- 
actions and the complex nature of hydrogen bonding. 

Models can be prepared either by introducing H in 
existing a-Si samples j8^,|9l|,|92| or by quenching a Si-H 
mixture from the melt |)3|j94| ]. The latter approach was 
employed by Tuttle and Adams (TA) f)3|| as well as Lan- 
zavecchia and Colombo (LC) [ p4[ . Both models contain 
216 Si atoms and either 24 or 26 H (i.e., about 10%). 
Using the SD-TB scheme, TA first quenched their cell 
from 1800 to 300 K at a rate of 10 15 K/s, equilibrated it 
for a short 0.5 ps, annealed it at 1200 K for 1.0 ps, then 
quenched it again to 300 K at a rate of 10 14 K/s. The 
final configuration (at 300 K) presents a relatively high 
concentration of defects: 1.5% of Si atoms are three-fold 
coordinated while 9.0% are five-fold coordinated, which 
leads to a gapless electronic density of states. For its 
part, hydrogen is fully bonded to Si, and found exclu- 
sively in monohydride configurations. The partial radial 
distribution functions, shown in Fig. |l3|, indicate that 
the silicon backbone is essentially identical to that of a 
pure a-Si network and that there is a complete absence 
of medium-range order in hydrogen. Most interestingly, 
the addition of hydrogen is sufficient to create local in- 
homogeneities in the cell. A small cavity, for example, is 
formed in the cell (Fig. |l4|), with H atoms "decorating" 
its internal surface. 

A similar calculation was carried out by LC using the 
GSP-TB model for Si and a similar one for Si-H and H-H 
interactions [B4f , as discussed in Section [II E 2 . The cool- 



ing rate in this case is four times smaller than that used 
by TA. Here again, monohydride complexes are found 
to be the most likely configuration for hydrogen. How- 
ever, LC also found some SiH2 complexes as well as small 
three- and four-atom H clusters. In spite of their large 
density of structural and electronic defects, these two 
models provide important hints on the structural modifi- 
cations of the amorphous phase which might be induced 
by H cavities, hydrogen chains and clusters. The dif- 
ference between TA and LC structures, in terms of the 
presence or absence of cavities and polyhydrides com- 
plexes, is not understood at the moment. It might be 
due to different modes of preparation or potentials, or to 
the fact that TA used a rescaled mass for H while LC 
employed the actual value. 

Starting from a "preformed" a-Si sample offers the ad- 
vantage that H can be introduced in a controlled manner 
and the release of electronic and elastic strain monitored, 
in particular through the removal of electronic traps. 
Gap states do not arise only from the presence of dangling 
bonds but also in strained environments, and hydrogen 
must therefore sometimes be forced in by breaking the 
network. For example, Fedders and Drabold |91| and 
Holender et al. [ |39[ simply removed overconstrained sili- 
con atoms and saturated the newly-formed dangly bonds 



with hydrogen. Evidently, such a procedure has little in 
common with the actual physical process, but neverthe- 
less proves useful insights on how hydrogen releases the 
strain and causes states to move deeper in the valence 
band. 

The Staebler-Wronski effect fl96|] , i.e., the degradation 
of the photoelectric conversion properties of the material 
under exposure to light, is the motivation behind much 
of the work on a-Si:H. This phenomenon is known to 
involve structural rearrangements following the absorp- 
tion of photons, and can in fact be reversed by annealing 
at sufficiently high temperature. It is however not clear 
if the electron emitted following the absorption of the 
photon plays an active role or if it merely transfers ki- 
netic energy to the network. A variety of approaches 
have been used, within the TBMD framework, to inves- 
tigate this problem. Using a 62-Si plus 5 or 7-H atom 
cell, Fedders |)7j has examined the effect, on energy and 
relaxation, of the state of charge of dangling bonds and 
finds the passivation energy to strongly depend on the lo- 
cal environment, fluctuating by as much as 1 eV in their 
limited sample. It also appears that charged dangling 
bonds would be more energetically favourable than neu- 
tral ones. It is difficult experimentally to establish the 
density of specific charge defects so that more numerical 
work is needed to confirm and refine the results of these 
calculations. 

A similar study was carried out by Biswas et al. 
(BLYB) H using the 60-atom (54 Si + 6 H) model of 
Guttman and Fong ]99|) . Here the cell was more com- 
pletely relaxed by equilibrating over several tens of ps at 
300 K while, in contrast, Fedders |)7j worked with con- 
figurations optimised locally at K. Adding or removing 
a H atom, BLYB followed the relaxation of the lattice. 
This relaxation takes place very rapidly, over a period of 
about 10 ps, and involves almost exclusively the Si atom 
bonded to the defect, with very similar behaviour for pos- 
itively and negatively charged defects. Park and Myles 
(PM) have proposed a "hot-spot" method for stimulat- 
ing relaxation |9^] . Starting with a K configuration, 
two atoms on a bond are given a burst of energy (2 eV 
in this case), and the dynamics is followed until the ex- 
citation reaches the boundary of the periodic box (300 
fs). Using the 60-atom Guttman- Fong model, Park and 
Myles have found, of all the bonds they tried, only one 
for which the hot-spot excitation leads to a new configu- 
ration, producing a dangling and a floating bonds, after 
jumping an activation barrier of about 2 eV. This is in 
agreement with BLYB's results where the localised and 
rapid relaxation is also an indication of a very stable con- 
figuration. In the context of the Staebler-Wronski effect, 
these results would support a very localised mechanism 
controlled by isolated bonds being broken and reforming. 
However, the stability of the lattice could also be due the 
its very small size; more simulations are needed in order 
to establish the microscopic origin of this effect. 

Other dynamical properties of a-Si:H have also been 
studied. Fedders and Drabold (FD), for instance, have 
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searched their model for conformational fluctuations by 
quenching snapshot configu ratio ns at regular intervals in 
time during a 600 K run 1 10C ] . They found that al- 
though the connectivity of the lattice is preserved, the 
quenched configurations are all slightly different, reveal- 
ing the presence of a large number of very small barri- 
ers between slightly different states. The most relevant 
barrier is that for the diffusion of H. This light atom is 
likely to move significantly in the mostly-empty struc- 
ture. However, as was the case for c-Si, the nature of the 
diffusion mechanism remains unclear. From simulations 
on a system containing 216 Si and 24 H atoms, Lanzavec- 
chia and Colombo found that H migrates through jumps 
between neighbouring dangling bonds via a metastablc 
Si-H-Si state . Si mila r observations were r epor ted by 
Tuttle and Adams [101] and Biswas et al. [102]: dur- 
ing diffusion, a hydrogen binds to an already four-fold- 
coordinated atom, giving rise, for a short period of time, 
to a floating bond; Biswas et al. estimate the formation 
energy of this floating bond to be in the range 1.3-2.3 
eV. 



C. Compounds 

Amorphous compounds, such as the III-V semiconduc- 
tors and the chalcogenides, are particularly challenging 
for simulations. First is the problem of constructing an 
appropriate set of interactions: not only are there more 
interaction types than in mono-atomic compounds, but 
also ionisation and charge transfer effects can become 
important, requiring special attention. Second, the po- 
tential energy surface is complicated by the introduction 
of a new dimension, namely the chemical identity of the 
constituent species. Thus, in a study of the liquid-to- 
amorphous transformation, for instance, the search for 
the ground state must seek to minimise both the topo- 
logical and the chemical disorder. As a consequence, the 
relaxation timescales available through TBMD simula- 
tions, which are already out of measure with experimen- 
tal timescales, become prohibitively long, except perhaps 
for some very ionic compounds, such as SiC>2, where even 
the liquid phase already exhibits almost perfect chemical 
order. 

In spite of this difficulty, most TBMD studies of amor- 
phous III-V materials proceeds via the usual melt-and- 
quench cycle. Indeed, there exists no "recipe" a la 
Wooten- Winer- Weaire to prepare amorphous compounds 
that are chemically ordered. Further, there exists no sat- 
isfactory classical potentials for these materials that can 
be used to bypass the computer-intensive TBMD melt- 
and-quench cycle. Nevertheless it is possible, through a 
combination of approaches, to generate very high-quality 
models, as will see below. 

TBMD can provide important insights on the physics 
of these materials in spite of the above limitations. 
This is demonstrated, for instance, in the recent simu- 



lation by Stumm and Drabold of GaN 102]. Using the 
Sankey-Drabold TBMD model (7§, Stumm and Drabold 
quenched a 64-atom cell of liquid GaN into the amor- 
phous phase at a rate of 1.3 x 10 15 K/s. Two densities 
were studied: one equal to that of wurtzite GaN and the 
other at 82% of this density. For both structures, the to- 
tal radial distribution function exhibits a deep minimum 
between the first and second-neighbour peaks despite a 
large concentration of three- fold atoms (37 and 66%, re- 
spectively). Contrary to what is seen in a-C, where a 
large density of three-fold atoms fills the electronic gap, 
there are no states which appear in the gap of either GaN 
unit cells. Amorphous GaN is thus predicted to be, in its 
own right, an interesting material for device applications. 

Studies of a-GaAs obtained by quenching from the liq- 
uid phase have been carried out by Molteni, Col omb o 
and Miglio (MCM) ||, and Seong and Lewis (SL) p5 |. 
This material has a relatively small ionic ity; the TB- 
potential propose d by Molteni et al. [ 104 1 , already dis- 
cussed in Section HID, deals with charge-transfer effects 



by introducing a screened electrostatic term, thus keep- 
ing the potential short range. In both cases, 64-atom cells 
were used and the cooling rates were similar (1.5 x 10 12 
K/s for MCM, 50% slower for SL). One important dif- 
ference however is in the choice of density: while MCM 
fixed it at the crystalline value, SL chose to optimise it; 
a density 3.2% smaller than the crystalline one was thus 
obtained (Fig. [l5]), in agreement with experiment — val- 
ues in the range 4.98-5.11 g/cm 3 have been reported; the 
crystalline value is 5.32 g/cm 3 . Optimisation of the den- 
sity also leads to significant, and perhaps unexpected, 
differences between the two structures: MCM find the 
bonding environments of As and Ga to be symmetric, 
that is the two species exhibit similar distributions of co- 
ordination defects; in contrast, SL observe sizable differ- 
ences between the two species, with Ga much more likely 
to be in a five-fold coordination state than As, which 
prefers to be in a three-fold state. Overall, however, the 
average coordination is almost precisely four — in fact 
slightly less, 3.94, certainly related to the fact that a- 
GaAs is less dense than c-GaAs. 

TBMD has also been used to generate m odels for the 
chalcogenide glasses. Cobb and co-workers [ 106, 107 1, for 
instance, have quenched 62- and 216-atom models of 
GeSe2 from the liquid phase using the Sankey-Drabold- 
TB scheme ]76| . The total simulation time for the 216- 
atom cell is 5 ps. In spite of this very short time, which 
leaves a large number of defects in the network, the sys- 
tem exhibits a well-defined optical gap of 1.72 eV, free 
of defect states; the static structure factor and vibra- 
tional density of states is also in reasonable agreement 
with experiment (Fig. |lj). Most (~85%) Ge atoms are 
fourfold coordinated with about 26% Ge having another 
Ge atom as one of its neighbours; a similar fraction of 
Se form homopolar bonds. Because of different coordi- 
nation, this leads to a total density of "wrong" bonds of 
10.8%. Cobb et al., furthermore, find that the localised 
electronic states are far from the band edges, thus leaving 
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a wide, state-free gap. 

Because GeSe2 is significantly less ionic than, e.g., 
SiC>2, the timescale needed to obtain a chemically- 
ordered GeSe2 glass is well beyond the reach of TBMD. 
This is true also of the III-V compound a-GaAs. As a 
result, the melt-and-quench simulations of a-GaAs men- 
tioned above lead to a density of defects which is large. 
Wrong bonds, for instance, are found to be in concentra- 
tions of 12.2 and 12.9% in the models of SL and MCM, 
respectively. While the exact number is not known from 
experiment, these values are most certainly on the high 
side. The full TBMD melt-and-quench cycle can be 
partly bypassed through a combination of approaches. 
As noted already, there exists not equivalent, in the 
case of III-V compounds, to the Wooten- Winer- Weaire 
model for a-Si; however, even though this model is "non- 
physical" (the topology is distorted in an ad hoc man- 
ner), it nevertheless yields structures in excellent agree- 
ment with experiment. In the same spirit — i.e., the 
end justifies the means — Mousseau and Lewis [ [73|j74| 
have proposed a scheme for a-GaAs based on empirical 
potentials that compares favourably with melt-quenched 
models, as we discuss now. 

Amorphous covalently-bonded semiconductors have 
traditionally been described in terms of CRN's. Polk's 
CRN |62), presumably appropriate to a-Si, consists of a 
collage of corner-sharing tetrahedra which preserves the 
ideal coordination of four. The Conncll-Tcmkin CRN 



[ 108 1 is similar, except for the additional constraint that 
odd-membered rings are not permitted. This model is 
likely relevant to chemically-ordered binary compounds 
since the presence of odd-membered rings necessarily im- 
ply the presence of wrong-bond defects. That a-GaAs is 
indeed akin to the Connell-Temkin model has been re- 
cently established by Mousseau and Lewis j73|,[74]]. This 
was done using a combination of approaches as follows: 
Starting with a 216-atom supercell in a random state, 
"zeroth-order" models were generated using Barkema 
and Mousseau's activation-relaxation technique (ART), 
the relaxation being carried out using a simple Stillinger- 
Weber-like potential with and without an additional term 
between like atoms; these structures were then relaxed 
using the MCM-TBMD model. Mousseau and Lewis 
have found the Connell-Temkin-type network to have sig- 
nificantly lower energy than the Polk-type network and 
therefore to provide a bet ter repre senta tion of a-GaAs 
as can be seen in Tables VII and VIIL However, the 



same structures with Ga and As replaced by Si and re- 
laxed using the GSP-TBMD showed no difference in en- 
ergy, indicating that both models have the same inherent 
strain level. This last result is rather surprising: the den- 
sity of odd-membered rings is apparently not determined 
by elastic constraints but, simply, by entropic considera- 
tions. 

Following the static minimisation stage, both GaAs 
models were relaxed at finite temperature (300 K for 7 
ps and 700 K for 8.8 ps) to assess their stability and 
compute their dynamical properties. It was found that 



the Polk CRN for GaAs is unstable and in fact distorts 
significantly at temperatures as low as 700 K while the 
Connell-Temkin CRN remains relatively unaffected by 
the annealing process. However, a detailed comparison 
of the structural, vibrational and electronic properties 
of the two networks reveals little qualitative differences 
between them, indicating that intermediate-range order 
plays a relatively minor role in determining the proper- 
ties of GaAs. It appears that only direct measurements 
of the density of wrong bonds can shed light on the ex- 
perimental nature of the structure of a-GaAs. 



D. Surfaces of amorphous semiconductors 

Surfaces clearly play an important role in the micro- 
electronic properties of devices and in particular influ- 
ence the growth process. However, there has been rela- 
tively few simulations of disordered-material surfaces (in- 
cluding growth) in the TB framework, possibly because 
these require, to start with, realistic (and large enough) 
sample of the bulk amorphous phase. It is important 
to note that the timescale for growth by atom deposi- 
tion is much shorter in computer models than in real 
experiments, so that it is difficult to extract quantitative 
information from such simulations. 



Kilian et al. [ 1 1 1 1 have used the ab initio TB scheme 
of Sankey and Drabold |76) to study the surface states 
of a-Si and a-Si:H. A 216-atom Wooten- Weiner- Weaire 
model of a-Si ||,||| was first relaxed using TBMD. The 
periodic boundary condition was then removed along the 
z direction and the bottom layer passivated with hydro- 
gen; the top layer, finally, was relaxed with and without 
hydrogen in order to study localised electronic surface 
states and the effect of hydrogenation on them. In order 
to do this, a local charge q(n, E) for atom n is defined, 
where E is the energy eigenstate. One may then calculate 
the mean square charge of a "layer" located at position 
z and containing N z atoms: 



q 2 (E,z)=J2l( 



Multiplying this quantity by N z gives a measure of the 
localisation of charge in the layer, Q 2 {E^ z); in this way, 
the decay of surface states with depth can be monitored. 
Figs. and |l| show some surface defects — not nec- 
essarily topological - and the corresponding spatial lo- 
calisation of charge. Hydrogenating the surface does not 
eliminate these states, however, contrary to what often 
happens in the bulk: Fig. [l8] shows that surface states 
are more localised by the addition of hydrogen but not 
otherwise affected. The study of other a-Si:H models 



[ 1 1 2 1 lead to similar conclusions: surface states seem to 



be difficult to eliminate. It should be noted however that, 
again because of computational limitations, the surfaces 
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may not be fully relaxed. It remains to be seen if de- 
tailed thermal annealing would decrease the influence of 
localised surface states. 



V. CONCLUDING REMARKS 

Tight-binding molecular dynamics has contributed sig- 
nificantly to our understanding of defects in semiconduc- 
tors. Because it is generally reliable, and yet economical 
from the computational viewpoint, the method consti- 
tutes an excellent complement to ab initio calculations. 
While the latter is limited to small unit cells, for which 
self-interaction of defects remains important, it is possi- 
ble with TBMD to deal with systems containing many 
hundreds of atoms in a fully-relaxed configuration. With 
the advent of O(N) methods, it becomes feasible to carry 
out TBMD studies of more complex defects involving 
many tens of atoms, a task which is at present beyond 
the reach of ab initio approaches. 

The study of amorphous semiconductors has also 
greatly benefited from TBMD calculations. With the 
wide range of local environments found in these materi- 
als, it is difficult to construct a fully satisfactory empirical 
potential; TB models have proved much more transfer- 
able. However, because of the computational cost of TB 
calculations, still much larger than classical potentials, 
and the difficulty of implementing O(N) methods for 
these materials, simulations have been limited to, princi- 
pally, 64-atom cells with some noteworthy exceptions up 
to 512-atom cells. Using a mixture of hands-on or empir- 
ical approaches and TBMD for creating the initial struc- 
tures can compensate for the increased computational 
cost. In spite of the biases that they possibly introduce, 
such mixed approaches appear to be the best, at present, 
for studying large, low-strain disordered structures. 

The examples presented in this article underline the 
important contributions of TBMD calculations to the 
study of semiconductors. As the method becomes more 
widespread and fundamentals are established, it can be 
expected that several of the questions left unanswered in 
this review will be addressed in the coming years. 
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TABLE I. Point defect formation energies (in eV) obtained using various TB models and comparison with the LDA results. 
The relaxation energies are given in parenthesis (when available). The LDA results for the two interstitials are not consistently 
relaxed because of the use of different computational schemes (e.g., supercell versus Green's functions). For the monovacancy 
and the divacancies, the LDA calculations allowed full relaxation. MC: Mercer and Chou jnj; KBWHS: Kwon et al. ]l2[ ; 
LKKVERYA: Lenosky et al. pf. 





GSP 


GSP GSP MC 


KBHS 


KBHS 


LKKVERYA LKKVERYA 


LDA 




TV = 64 


TV = 216 TV = 512 TV = 64 


TV = 64 


TV = 216 


TV = 64 


TV = 216 


TV = 64 


References 


[9,17] 


[i | [hi 




[12| 


[13| 


[ll 


14 ie] 


Defect 





Monoacancy 
T interstitial 
H interstitial 
Divacancy 
Split divacancy 
Frenkel pair 



3.67 (1.35) 3.96 (1.86) 4.12 (1.91) 
4.39 (2.13) 4.40 (2.14) 4.41 (2.14) 
5.78 (4.30) 5.90 (4.42) 5.93 (4.42) 

5.68 (1.73) 

6.54 (2.43) 

6.55 (3.89) 



3.76 3.46 (1.26) 3.93 (1.64) 3.40 (0.24) 
4.95 3.61 (0.51) 4.42 (0.49) 3.55 (0.59) 
4.75 (1.17) 5.13 (1.23) 3.56 (0.89) 



3.78 (0.20) 



3.3-4.3 (0.4-0.6) 
3.7-4.8 (0.1-0.2) 
4.3-5.0 (0.6-1.1) 
4.32 (0.27) 
5.90 (0.57) 



TABLE II. Energy levels of various neutral defects in their relaxed configurations, in eV, measured with respect to the 
valence-band maximum. GSP-TBMD: Results from Ref. 0; LDA: Results from Ref. jl6). 



Defect 


GSP-TBMD 


LDA 


Monoacancy 


0.76 


0.23 


T interstitial 


0.52, 0.12 




Divacancy 


1.00, 0.46 


0.04 


Split divacancy 


0.76, 0.46, 0.34 


0.08 


Frenkel pair 


0.76, 0.52 





TABLE III. s-like (A\ symmetry) deep energy levels in GaP and Si, measured with respect to the top of the valence band; 
C.B. indicates levels lying in the conduction band. TB: TB model of Hjalmarson et al. (So), which does not include relaxation. 
TBMD: TB model of Hjalmarson et al. [B(J including relaxation, after Li and Myles, Ref. EjJ, where experimental references 



can also be found. The equilibrium bond lengths 


are 2.36 and 2.35 A for GaP and Si, respectively. 




System 


Bond length 


TB 


TBMD 


Experiment 




A 


eV 


eV 


eV 


GaP:N 


2.15 


2.10 


2.25 


2.34 


GaP:0 


2.78 


1.85 


1.70 


1.46 


GaP:P Ga 


2.28 


1.03 


1.09 


1.10 


GaP:Ge 


2.24 


1.85 


1.95 


2.16 


GaP:Se 


2.22 


2.32 


C.B. 




Si:S 


2.23 


0.58 


0.63 


0.85 


Si:Se 


2.03 


0.65 


0.83 


0.86 


Si:Te 


2.60 


1.12 


1.05 


1.01 


Si:C 


1.98 


1.09 


C.B. 
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TABLE IV. Nearest-neighbour breathing-mode displacements and local volume changes for vacancies and antisites in GaAs, 
after Seong and Lewis [Q. The average relaxations are given both in A and relative to the bond distance of bulk GaAs (2.45 
A). + and — refer to outward and inward relaxation, respectively. AV = V — Vo is the change in volume of the defect (i.e., the 
volume of the tetrahedron formed by the four nearest neighbours) resulting from relaxation. 



Defect 


1 


Breathing 
2 


(A) 

3 


4 


Average (A) 


Average (%) 


AV/V (%) 


vL 


-0.31 


-0.35 


-0.31 


-0.31 


-0.32 


-13.1 


-34.3 




-0.33 


-0.33 


-0.31 


-0.30 


-0.32 


-13.0 


-34.1 




-0.32 


-0.32 


-0.32 


-0.29 


-0.32 


-12.9 


-33.9 




-0.32 


-0.32 


-0.32 


-0.30 


-0.31 


-12.8 


-33.7 


vl~ 


-0.31 


-0.31 


-0.31 


-0.31 


-0.31 


-12.8 


-33.7 


vL 


0.38 


0.38 


0.38 


0.38 


0.38 


15.5 


54.0 


vL 


0.44 


0.44 


0.44 


-0.49 


0.21 


8.5 


24.2 




0.39 


0.39 


0.39 


-0.50 


0.17 


6.8 


18.9 


vl; 


0.45 


-0.42 


0.45 


-0.42 


0.02 


0.8 


-3.3 


Ga As 


0.49 


0.49 


0.21 


0.21 


0.35 


14.4 


50.0 


Ga As 


0.48 


0.17 


0.13 


0.13 


0.23 


9.2 


30.0 


G a As 


0.48 


0.07 


0.07 


0.07 


0.17 


7.0 


22.2 


Ga As 


0.11 


0.02 


0.02 


0.02 


0.04 


1.8 


5.5 


Ga As 


-0.01 


-0.01 


-0.01 


-0.01 


-0.01 


-0.3 


-1.0 


A4t 
As£ a 


0.06 


0.06 


0.06 


0.06 


0.06 


2.5 


7.6 


0.06 


0.06 


0.06 


0.06 


0.06 


2.4 


7.5 


As& a 


0.06 


0.06 


0.06 


0.06 


0.06 


2.4 


7.5 


As Ga 


0.06 


0.06 


0.06 


0.06 


0.06 


2.5 


7.8 



TABLE V. Nearest-neighbour pairing-mode displacements for As vacancies and Ga antisites, after Seong and Lewis 



Defect 




Pairing 


1(A) 






Pairing 2 (A) 






1 


2 


3 


4 


1 


2 


3 


4 




























vl 


-0.01 


-0.01 


0.03 





0.02 


0.02 










-0.01 


-0.01 


0.02 

















vl: 


0.01 


-0.01 


0.01 


-0.01 


-0.01 


0.01 


-0.01 


0.01 


Ga As 


-0.02 


-0.02 


-0.21 


-0.21 














Ga As 


0.03 


-0.01 


-0.12 


-0.12 








-0.03 


0.03 


Ga As 





0.08 


-0.04 


-0.04 








-0.06 


0.06 


Ga As 





0.05 


-0.03 


-0.03 








-0.05 


0.05 


Ga As 



























TABLE VI. Coordination and bond angle distribution for a variety of o-Si models. CP: Stich et al. tab initio quench) |p9j; 
ML: Mousseau and Lewis (empirical construction, TB relaxation) frs} ]; KL: Kim and Lee (TB quench) |37j; DFSD: Drabold et 
al. (TB quench) Q; SC: Servalli and Colombo (TB quench) [f|. 





CP 


ML 


KL 


DFSD 


SCA 


SCB 


Z = 3 


0.002 


0.032 


0.032 


0.063 


0.049 


0.009 


Z = 4 


0.966 


0.954 


0.828 


0.877 


0.904 


0.958 


Z = 5 


0.032 


0.014 


0.125 


0.060 


0.047 


0.033 


Z = 6 








0.016 









<Z> 


4.03 


3.98 


4.28 


4.00 


3.99 


4.02 


ro 


3.0 


3.0 


2.75 


2.70 


2.82 


2.82 


6 


106.7 


109.4 


108.3 








A6 


16.3 


9.4 


15.5 








# atoms 


64 


216 


64 


63 


64 


64 
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TABLE VII. Structural characteristics of various models for a-GaAs. Distribution of coordination numbers, Z (and near- 
est-neighbour cutoff distance, rjyjv), first nearest-neighbour distance, r\ (and width, Ar^), density of wrong bonds, and width 
of the bond-angle distribution, A8. CRN-P: Polk-type continuous random network of Mousseau and Lewis j7^, CRN-CT: 
Connell-Temkin-type continuous random network of Mousseau and Le wis |73| , |74) ; SL — TB simulations of Ref. [ 105[ | ; MCM — 
TB simulations of Ref. j35|; CP — Car-Parrinello simulations of Ref. | 109[ 







CRN-P 




CRN-CT 


SL 


MCM 


CP 




K 


300 K 


K 


300 K 


K 


K 


K 


Z = 3 


0.046 


0.128 


0.051 


0.118 


0.242 


0.14 


0.219 


Z = 4 


0.954 


0.845 


0.944 


0.830 


0.598 


0.66 


0.781 


Z = 5 





0.026 


0.005 


0.045 


0.129 


0.18 





Z = 6 





0.001 





0.004 


0.024 







Z= 7 





0.000 





0.002 


0.007 







<Z> 


3.95 


3.90 


3.95 


3.95 


3.94 


4.09 


3.83 


r NN (A) 


3.0 


3.1 


3.0 


3.1 


3.0 


3.0 


2.8 


n (A) 


2.508 


2.505 


2.517 


2.507 








An (A) 


0.075 


0.117 


0.073 


0.103 








Wrong bonds (%) 


14.1 


14.2 


3.9 


5.2 


12.2 


12.9 


10.0 


A 6* (deg.) 


11.0 


14.1 


10.8 


15.0 


17.0 


17.0 





TABLE VIII. Energy (eV/atom) of the Polk-type (CRN-P) and Connell-Temkin-type (CRN-CT) 
a-GaAs continuous random networks of Mousseau and Lewis, relaxed with TB potentials at K [|T3|J 
give the results from the TB-MD simulations of Seong and Lewis (SL), Ref. 105 1. 



models for a-Si and 
. For GaAs, we also 



Network 


Si 


GaAs 


GaAs (SL) 


CRN-P 


-13.172 


-13.450 




CRN-CT 


-13.163 


-13.561 


-13.450 


Crystal 


-13.389 


-13.802 


-13.802 
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FIG. 1. Density of states near the band gap for various defects in crystalline Si. (From Ref. |17 
permission) . 



reproduced by kind 
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FIG. 2. Number of dangling bonds N&b (top panel), formation energy Ef (middle panel) and binding energy Eb (bottom 
panel) for clusters of vacancies in Si as a function of size n. The full lines and diamond symbols are for HRC clusters, while 
dashed lines and plus symbols are for SPC clusters. (From Ref. [p7[; reproduced by kind permission). 
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(a) Ga A! 



o 



(b) Ga A 





°6° 



FIG. 3. Local relaxed atomic configuration for (a) the negative and (b) the neutral Ga antisite in GaAs. (From Ref. 
reproduced by kind permission). 
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FIG. 5. Defect formation energies in GaAs vs A(iQ a for three different values of the electron chemical potential: valence-band 
maximum (VBM), midgap (E g /2), and conduction band minimum (CBM). The lower and upper limits of the A^iGa— range 
correspond to As and Ga-rich regions, respectively. (From Ref. j34|; reproduced by kind permission). 
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FIG. 6. Arrhenius plot of the diffusion constant for vacancies and interstitials in crystalline silicon. (From Ref. 
reproduced by kind permission). 
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FIG. 8. Airhenius plot of the diffusion constant of H in c-Si. The open squares are the TBMD results of Panzarini and 
Colombo H] and the open diamonds are the ab initio results of Buda et al. |55) ; also given are the corresponding experimental 
data (cf. [bij for references; reproduced by kind permission). 
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FIG. 9. Diffusion path of H in c-Si. (From Ref. [Q; reproduced by kind permission). 
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FIG. 10. TB electronic density ol states (spin/eV/atom) for various a-Si models obtained using empirical potentials, 
(a) 216-atom Wooten-Weiner-Weaire (WWW) model (no coordination defects); (b) 216-atom WWW model relaxed with a 
Stillinger- Weber potential (no coordination defects); (c) 216-atom WWW/MD model (two three-fold defects); (d) 1728-atom 
model (4% five-fold and 2% three-fold defects); (e) 1728-atom WWW/MD model (no coordination defects); (f) 13824-atom 
WWW/MD model (4.5% five-fold and 2% three-fold defects). (From Ref. [^o|; reproduced by kind permission). 
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FIG. 11. Electronic states in the band-gap region of a-Si: (a) Electronic density of states; inset: total density of states, (b) 



Inverse participation ratio (a measure of localisation) vs ener 



of doping (computed in the Kubo formalism). (From Ref. |81 , reproduced by kind permission.) 



5y for states close to the gap. (c) DC conductivity as a function 
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FIG. 13. Partial radial distribution functions for hydrogenated amorphous silicon. Full line: experimental measurements of 
Menelle [^H); dashed line: 216 Si + 26 H model of Tuttle and Adams [Q; dotted line: ab initio simulation of the (61 Si + 11 
H)-atom cell of Buda et al. 110]. (From Ref. tea]; reproduced by kind permission). 
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FIG. 14. Ellipsoidal cavity in a 242-atom a-Si:H model containing 26 H, prepared with MD from a liquid configuration. 
Eight H atoms bonded to monohydride Si form the surface of this cavity of dimensions 10 x 5 x 5 A 3 . (From Ref. ji)^]; reproduced 
by kind permission). 
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FIG, 15. Total energy vs. lattice constant at K for crystalline GaAs and the 64-atom amorphous samples of Seong and 
Lewis 1 105 1 . (Reproduced by kind permission). 
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FIG. 16. Vibrational density ol states fo r am orphous GeS e2. T he solid curve is obtained from the 216-atom cell of Cobb et 
al. |l07| and the dots are from experiment 106 1. (From Ref. 107 1; reproduced by kind permission). 
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FIG. 18. Spatial charge localisation as a function of distance from the bottom layer for defects (b) of Fig. [n], before (top) 
and after (bottom) hydrogenation of the surface. It can be seen t hat t he electronic defects do not disappear upon hydrogenation 
but merely localise somewhat closer to the surface. (From Ref. [Ill ; reproduced by kind permission). 
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